home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dggbal.f < prev    next >
Text File  |  1996-11-04  |  13KB  |  462 lines
  1.       SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
  2.      $                   RSCALE, WORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOB
  11.       INTEGER            IHI, ILO, INFO, LDA, LDB, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), LSCALE( * ),
  15.      $                   RSCALE( * ), WORK( * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DGGBAL balances a pair of general real matrices (A,B).  This
  22. *  involves, first, permuting A and B by similarity transformations to
  23. *  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
  24. *  elements on the diagonal; and second, applying a diagonal similarity
  25. *  transformation to rows and columns ILO to IHI to make the rows
  26. *  and columns as close in norm as possible. Both steps are optional.
  27. *
  28. *  Balancing may reduce the 1-norm of the matrices, and improve the
  29. *  accuracy of the computed eigenvalues and/or eigenvectors in the
  30. *  generalized eigenvalue problem A*x = lambda*B*x.
  31. *
  32. *  Arguments
  33. *  =========
  34. *
  35. *  JOB     (input) CHARACTER*1
  36. *          Specifies the operations to be performed on A and B:
  37. *          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
  38. *                  and RSCALE(I) = 1.0 for i = 1,...,N.
  39. *          = 'P':  permute only;
  40. *          = 'S':  scale only;
  41. *          = 'B':  both permute and scale.
  42. *
  43. *  N       (input) INTEGER
  44. *          The order of the matrices A and B.  N >= 0.
  45. *
  46. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  47. *          On entry, the input matrix A.
  48. *          On exit,  A is overwritten by the balanced matrix.
  49. *          If JOB = 'N', A is not referenced.
  50. *
  51. *  LDA     (input) INTEGER
  52. *          The leading dimension of the array A. LDA >= max(1,N).
  53. *
  54. *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
  55. *          On entry, the input matrix B.
  56. *          On exit,  B is overwritten by the balanced matrix.
  57. *          If JOB = 'N', B is not referenced.
  58. *
  59. *  LDB     (input) INTEGER
  60. *          The leading dimension of the array B. LDB >= max(1,N).
  61. *
  62. *  ILO     (output) INTEGER
  63. *  IHI     (output) INTEGER
  64. *          ILO and IHI are set to integers such that on exit
  65. *          A(i,j) = 0 and B(i,j) = 0 if i > j and
  66. *          j = 1,...,ILO-1 or i = IHI+1,...,N.
  67. *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  68. *
  69. *  LSCALE  (output) DOUBLE PRECISION array, dimension (N)
  70. *          Details of the permutations and scaling factors applied
  71. *          to the left side of A and B.  If P(j) is the index of the
  72. *          row interchanged with row j, and D(j)
  73. *          is the scaling factor applied to row j, then
  74. *            LSCALE(j) = P(j)    for J = 1,...,ILO-1
  75. *                      = D(j)    for J = ILO,...,IHI
  76. *                      = P(j)    for J = IHI+1,...,N.
  77. *          The order in which the interchanges are made is N to IHI+1,
  78. *          then 1 to ILO-1.
  79. *
  80. *  RSCALE  (output) DOUBLE PRECISION array, dimension (N)
  81. *          Details of the permutations and scaling factors applied
  82. *          to the right side of A and B.  If P(j) is the index of the
  83. *          column interchanged with column j, and D(j)
  84. *          is the scaling factor applied to column j, then
  85. *            LSCALE(j) = P(j)    for J = 1,...,ILO-1
  86. *                      = D(j)    for J = ILO,...,IHI
  87. *                      = P(j)    for J = IHI+1,...,N.
  88. *          The order in which the interchanges are made is N to IHI+1,
  89. *          then 1 to ILO-1.
  90. *
  91. *  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)
  92. *
  93. *  INFO    (output) INTEGER
  94. *          = 0:  successful exit
  95. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  96. *
  97. *  Further Details
  98. *  ===============
  99. *
  100. *  See R.C. WARD, Balancing the generalized eigenvalue problem,
  101. *                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
  102. *
  103. *  =====================================================================
  104. *
  105. *     .. Parameters ..
  106.       DOUBLE PRECISION   ZERO, HALF, ONE
  107.       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
  108.       DOUBLE PRECISION   THREE, SCLFAC
  109.       PARAMETER          ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
  110. *     ..
  111. *     .. Local Scalars ..
  112.       INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
  113.      $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
  114.      $                   M, NR, NRP2
  115.       DOUBLE PRECISION   ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
  116.      $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
  117.      $                   SFMIN, SUM, T, TA, TB, TC
  118. *     ..
  119. *     .. External Functions ..
  120.       LOGICAL            LSAME
  121.       INTEGER            IDAMAX
  122.       DOUBLE PRECISION   DDOT, DLAMCH
  123.       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
  124. *     ..
  125. *     .. External Subroutines ..
  126.       EXTERNAL           DAXPY, DSCAL, DSWAP, XERBLA
  127. *     ..
  128. *     .. Intrinsic Functions ..
  129.       INTRINSIC          ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
  130. *     ..
  131. *     .. Executable Statements ..
  132. *
  133. *     Test the input parameters
  134. *
  135.       INFO = 0
  136.       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
  137.      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
  138.          INFO = -1
  139.       ELSE IF( N.LT.0 ) THEN
  140.          INFO = -2
  141.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  142.          INFO = -4
  143.       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  144.          INFO = -5
  145.       END IF
  146.       IF( INFO.NE.0 ) THEN
  147.          CALL XERBLA( 'DGGBAL', -INFO )
  148.          RETURN
  149.       END IF
  150. *
  151.       K = 1
  152.       L = N
  153. *
  154. *     Quick return if possible
  155. *
  156.       IF( N.EQ.0 )
  157.      $   RETURN
  158. *
  159.       IF( LSAME( JOB, 'N' ) ) THEN
  160.          ILO = 1
  161.          IHI = N
  162.          DO 10 I = 1, N
  163.             LSCALE( I ) = ONE
  164.             RSCALE( I ) = ONE
  165.    10    CONTINUE
  166.          RETURN
  167.       END IF
  168. *
  169.       IF( K.EQ.L ) THEN
  170.          ILO = 1
  171.          IHI = 1
  172.          LSCALE( 1 ) = ONE
  173.          RSCALE( 1 ) = ONE
  174.          RETURN
  175.       END IF
  176. *
  177.       IF( LSAME( JOB, 'S' ) )
  178.      $   GO TO 190
  179. *
  180.       GO TO 30
  181. *
  182. *     Permute the matrices A and B to isolate the eigenvalues.
  183. *
  184. *     Find row with one nonzero in columns 1 through L
  185. *
  186.    20 CONTINUE
  187.       L = LM1
  188.       IF( L.NE.1 )
  189.      $   GO TO 30
  190. *
  191.       RSCALE( 1 ) = 1
  192.       LSCALE( 1 ) = 1
  193.       GO TO 190
  194. *
  195.    30 CONTINUE
  196.       LM1 = L - 1
  197.       DO 80 I = L, 1, -1
  198.          DO 40 J = 1, LM1
  199.             JP1 = J + 1
  200.             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
  201.      $         GO TO 50
  202.    40    CONTINUE
  203.          J = L
  204.          GO TO 70
  205. *
  206.    50    CONTINUE
  207.          DO 60 J = JP1, L
  208.             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
  209.      $         GO TO 80
  210.    60    CONTINUE
  211.          J = JP1 - 1
  212. *
  213.    70    CONTINUE
  214.          M = L
  215.          IFLOW = 1
  216.          GO TO 160
  217.    80 CONTINUE
  218.       GO TO 100
  219. *
  220. *     Find column with one nonzero in rows K through N
  221. *
  222.    90 CONTINUE
  223.       K = K + 1
  224. *
  225.   100 CONTINUE
  226.       DO 150 J = K, L
  227.          DO 110 I = K, LM1
  228.             IP1 = I + 1
  229.             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
  230.      $         GO TO 120
  231.   110    CONTINUE
  232.          I = L
  233.          GO TO 140
  234.   120    CONTINUE
  235.          DO 130 I = IP1, L
  236.             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
  237.      $         GO TO 150
  238.   130    CONTINUE
  239.          I = IP1 - 1
  240.   140    CONTINUE
  241.          M = K
  242.          IFLOW = 2
  243.          GO TO 160
  244.   150 CONTINUE
  245.       GO TO 190
  246. *
  247. *     Permute rows M and I
  248. *
  249.   160 CONTINUE
  250.       LSCALE( M ) = I
  251.       IF( I.EQ.M )
  252.      $   GO TO 170
  253.       CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
  254.       CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
  255. *
  256. *     Permute columns M and J
  257. *
  258.   170 CONTINUE
  259.       RSCALE( M ) = J
  260.       IF( J.EQ.M )
  261.      $   GO TO 180
  262.       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
  263.       CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
  264. *
  265.   180 CONTINUE
  266.       GO TO ( 20, 90 )IFLOW
  267. *
  268.   190 CONTINUE
  269.       ILO = K
  270.       IHI = L
  271. *
  272.       IF( ILO.EQ.IHI )
  273.      $   RETURN
  274. *
  275.       IF( LSAME( JOB, 'P' ) )
  276.      $   RETURN
  277. *
  278. *     Balance the submatrix in rows ILO to IHI.
  279. *
  280.       NR = IHI - ILO + 1
  281.       DO 200 I = ILO, IHI
  282.          RSCALE( I ) = ZERO
  283.          LSCALE( I ) = ZERO
  284. *
  285.          WORK( I ) = ZERO
  286.          WORK( I+N ) = ZERO
  287.          WORK( I+2*N ) = ZERO
  288.          WORK( I+3*N ) = ZERO
  289.          WORK( I+4*N ) = ZERO
  290.          WORK( I+5*N ) = ZERO
  291.   200 CONTINUE
  292. *
  293. *     Compute right side vector in resulting linear equations
  294. *
  295.       BASL = LOG10( SCLFAC )
  296.       DO 240 I = ILO, IHI
  297.          DO 230 J = ILO, IHI
  298.             TB = B( I, J )
  299.             TA = A( I, J )
  300.             IF( TA.EQ.ZERO )
  301.      $         GO TO 210
  302.             TA = LOG10( ABS( TA ) ) / BASL
  303.   210       CONTINUE
  304.             IF( TB.EQ.ZERO )
  305.      $         GO TO 220
  306.             TB = LOG10( ABS( TB ) ) / BASL
  307.   220       CONTINUE
  308.             WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
  309.             WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
  310.   230    CONTINUE
  311.   240 CONTINUE
  312. *
  313.       COEF = ONE / DBLE( 2*NR )
  314.       COEF2 = COEF*COEF
  315.       COEF5 = HALF*COEF2
  316.       NRP2 = NR + 2
  317.       BETA = ZERO
  318.       IT = 1
  319. *
  320. *     Start generalized conjugate gradient iteration
  321. *
  322.   250 CONTINUE
  323. *
  324.       GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
  325.      $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
  326. *
  327.       EW = ZERO
  328.       EWC = ZERO
  329.       DO 260 I = ILO, IHI
  330.          EW = EW + WORK( I+4*N )
  331.          EWC = EWC + WORK( I+5*N )
  332.   260 CONTINUE
  333. *
  334.       GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
  335.       IF( GAMMA.EQ.ZERO )
  336.      $   GO TO 350
  337.       IF( IT.NE.1 )
  338.      $   BETA = GAMMA / PGAMMA
  339.       T = COEF5*( EWC-THREE*EW )
  340.       TC = COEF5*( EW-THREE*EWC )
  341. *
  342.       CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
  343.       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
  344. *
  345.       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
  346.       CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
  347. *
  348.       DO 270 I = ILO, IHI
  349.          WORK( I ) = WORK( I ) + TC
  350.          WORK( I+N ) = WORK( I+N ) + T
  351.   270 CONTINUE
  352. *
  353. *     Apply matrix to vector
  354. *
  355.       DO 300 I = ILO, IHI
  356.          KOUNT = 0
  357.          SUM = ZERO
  358.          DO 290 J = ILO, IHI
  359.             IF( A( I, J ).EQ.ZERO )
  360.      $         GO TO 280
  361.             KOUNT = KOUNT + 1
  362.             SUM = SUM + WORK( J )
  363.   280       CONTINUE
  364.             IF( B( I, J ).EQ.ZERO )
  365.      $         GO TO 290
  366.             KOUNT = KOUNT + 1
  367.             SUM = SUM + WORK( J )
  368.   290    CONTINUE
  369.          WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
  370.   300 CONTINUE
  371. *
  372.       DO 330 J = ILO, IHI
  373.          KOUNT = 0
  374.          SUM = ZERO
  375.          DO 320 I = ILO, IHI
  376.             IF( A( I, J ).EQ.ZERO )
  377.      $         GO TO 310
  378.             KOUNT = KOUNT + 1
  379.             SUM = SUM + WORK( I+N )
  380.   310       CONTINUE
  381.             IF( B( I, J ).EQ.ZERO )
  382.      $         GO TO 320
  383.             KOUNT = KOUNT + 1
  384.             SUM = SUM + WORK( I+N )
  385.   320    CONTINUE
  386.          WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
  387.   330 CONTINUE
  388. *
  389.       SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
  390.      $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
  391.       ALPHA = GAMMA / SUM
  392. *
  393. *     Determine correction to current iteration
  394. *
  395.       CMAX = ZERO
  396.       DO 340 I = ILO, IHI
  397.          COR = ALPHA*WORK( I+N )
  398.          IF( ABS( COR ).GT.CMAX )
  399.      $      CMAX = ABS( COR )
  400.          LSCALE( I ) = LSCALE( I ) + COR
  401.          COR = ALPHA*WORK( I )
  402.          IF( ABS( COR ).GT.CMAX )
  403.      $      CMAX = ABS( COR )
  404.          RSCALE( I ) = RSCALE( I ) + COR
  405.   340 CONTINUE
  406.       IF( CMAX.LT.HALF )
  407.      $   GO TO 350
  408. *
  409.       CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
  410.       CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
  411. *
  412.       PGAMMA = GAMMA
  413.       IT = IT + 1
  414.       IF( IT.LE.NRP2 )
  415.      $   GO TO 250
  416. *
  417. *     End generalized conjugate gradient iteration
  418. *
  419.   350 CONTINUE
  420.       SFMIN = DLAMCH( 'S' )
  421.       SFMAX = ONE / SFMIN
  422.       LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
  423.       LSFMAX = INT( LOG10( SFMAX ) / BASL )
  424.       DO 360 I = ILO, IHI
  425.          IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
  426.          RAB = ABS( A( I, IRAB+ILO-1 ) )
  427.          IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDA )
  428.          RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
  429.          LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
  430.          IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
  431.          IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
  432.          LSCALE( I ) = SCLFAC**IR
  433.          ICAB = IDAMAX( IHI, A( 1, I ), 1 )
  434.          CAB = ABS( A( ICAB, I ) )
  435.          ICAB = IDAMAX( IHI, B( 1, I ), 1 )
  436.          CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
  437.          LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
  438.          JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
  439.          JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
  440.          RSCALE( I ) = SCLFAC**JC
  441.   360 CONTINUE
  442. *
  443. *     Row scaling of matrices A and B
  444. *
  445.       DO 370 I = ILO, IHI
  446.          CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
  447.          CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
  448.   370 CONTINUE
  449. *
  450. *     Column scaling of matrices A and B
  451. *
  452.       DO 380 J = ILO, IHI
  453.          CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
  454.          CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
  455.   380 CONTINUE
  456. *
  457.       RETURN
  458. *
  459. *     End of DGGBAL
  460. *
  461.       END
  462.